A vortex dipole in a trapped two-dimensional Bose-Einstein condensate 
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We study the conservative dynamics and stationary configurations of a vortex-antivortex pair 
in a harmonically trapped two-dimensional Bose-Einstein condensate. We establish the conceptual 
framework for understanding the stationary states and the topological defect trajectories, through 
considerations of different mechanisms of vortex motion and the bifurcation of soliton-like stationary 
solutions. Our insights are based on Lagrangian-based variational calculations, numerical solutions 
of both the time-dependent and time-independent Gross-Pitaevskii equations, and exact solutions 
for the non-interacting case. 



I. INTRODUCTION 

Vortex-antivortex pairs or vortex dipoles are ubiqui- 
tous in two-dimensional (2D) Bose condensates. They 
are the most natural excitations when a 2D fluid flows 
past a barrier or through a disordered medium. Such 
creation of vortices and antivortices has been seen ex- 
plicitly in recent experiments with trapped atomic con- 
densates [H. In 2D, vortex dipoles can be thermally cre- 
ated. This is important for the thermodynamics of a 2D 
superfluid, e.g., for the Kosterlitz-Thouless phase tran- 
sition recently studied experimentally |2j. Given their 
pervasiveness and the current experimental interest with 
harmonically trapped condensates, a thorough study of 
vortex dipoles in trapped condensates is clearly impor- 
tant. 

In this Article, we use a combination of methods to 
analyze vortex-antivortex states and dynamics in a cir- 
cularly trapped 2D condensate, as described by the time- 
dependent Gross-Pitaevskii equation (GPE). First, we 
use a variational approach to the GPE that can be solved 
almost analytically, providing qualitative insight. These 
insights are explored and complemented by explicit nu- 
merical evolutions of the time-dependent GPE, exact 
solutions for the non-interacting case, and a numerical 
procedure for accessing stationary solutions of the time- 
independent GPE. Our results indicate that the dynam- 
ics of a vortex-antivortex pair in a trap is unexpectedly 
rich. In this Article we focus on a subset of possible dy- 
namics, namely, the pair trajectory for cases where the 
pair is initially placed symmetrically on opposite sides of 
the trap center. 

The physics behind our characteristic pair trajectories 
involves two effects. First, a vortex dipole in a large fluid 
body is a propagating object, i.e., the vortex and an- 
tivortex propel each other in a direction perpendicular 
to the line joining them. This has been studied in con- 
densates 0, [ij and is also common experience for ordi- 
nary fluids. In a trapped non-rotating condensate, there 
is an additional effect. An off-center vortex is known 
to precess around the center of the trap 0, Q; we can 
regard this as inhomogeneity-driven motion. For a vor- 
tex dipole in a trap, each of the pair is driven by the 
inhomogeneity, in a direction opposite to the mutually- 



driven motion. One of the two types of motion, mutually 
driven or inhomogeneity-driven, can dominate, depend- 
ing on the vortex-antivortex distance. This leads to the 
remarkable fact that a vortex dipole with the same sense 
(same dipole direction) can be propagating in one of two 
opposite directions, depending on the dipole separation. 
A similar situation holds for vortex rings in 3D trapped 
condensates @, HI ■ Another result of the competition de- 
scribed above is that the two effects can in some cases 
exactly balance each other, leading to stationary config- 
urations of vortex pairs. In this Article, we will examine 
in detail effects of this competition, for a vortex dipole 
in a circularly trapped condensate. 

The conservative dynamics presented here should be 
accessible experimentally, e.g., by creating a vortex- 
antivortex pair by phase imprinting, analogous to the 
first creation of a vortex in a laser-cooled condensate 0] . 
In contrast, the recent experiments P, deal with dissi- 
pative physics involving creation and destruction of vor- 
tex dipoles. While we have omitted dissipative consid- 
erations in this study, we will see that the conservative 
dynamics is extremely rich by itself. 

Theoretically, vortex dipoles in trapped condensates 
have appeared as a peripheral topic in numerical studies 
of 3D vortex rings 0, H, [HI [Hi Q3 , which may be re- 
garded as 3D analogs of 2D vortex dipoles. Refs. (l3l.[l4j| 
have considered stationary vortex dipole states, which 
we will further examine and interpret, Ref. [l5| has stud- 
ied energetics, while Ref. [16 | pr esented some density dy- 
namics. Very recently, Ref. [l7| has studied some vortex 
dipole dynamics in the weak-interaction region. 

Sec. HI presents results from our variational calcula- 
tion. Although the variational results cannot be expected 
to be quantitatively accurate, the calculation provides 
physical intuition. In particular, the most prominent 
feature of the vortex dipole trajectories, involving each 
singularity revolving around a stationary point, already 
emerges from this calculation. In Sec. IHII we present rel- 
evant exact results for the non-interacting case, focusing 
again on defect trajectories. Sec. IIVI presents numerical 
solutions of the time-independent GPE. This provides 
valuable insight into the differences observed in the dy- 
namics for small interactions and large interactions. In 
Sec. |V] we present direct numerical simulations of the 
time-dependent GPE. The earlier calculations (Sees. [Til 
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FIG. 1: (Color online.) Initial pair velocities versus dipole 
separation. Thick-solid: Gaussian ansatz, g — 0.1. Dashed: 
a single vortex, same ansatz, same g — 0.1. Dotted line: non- 
interacting condensate, ^xj — 2/xa- Inset: Initial velocities 
versus Xd, Solid, dashed and dotted-dashed lines are from 
Thomas-Fermi ansatz, with g = 1, 10, and 100. Dotted line 
is from gaussian ansatz, g — 1. All quantities are plotted in 
trap units. 

IIII1 IIVI) allow us to interpret the results of the time- 
dependent simulations; Sees. I VII and I VIII present the in- 
terpretations and point out open questions. 

Gross- Pitaevskii Equation. We study a low- 

temperature Bose-Einstein condensate in a harmonic 
trap. From the outset we start with a two-dimensional 
system confined to the xy plane, assuming that the dy- 
namics in the z direction is frozen out by tight confine- 
ment which we do not treat explicitly. We will consider 
circular traps, and measure lengths in units of trap os- 
cillator length and time in units of inverse trapping fre- 
quency. The condensate dynamics is given by the time- 
dependent Gross-Pitaevskii equation [3, [l9[ : 

*^ = -5V 2 ^ + H r (r)^ + ff |V'rV. (i) 

The trap potential is Vtr = \{x 2 + y 2 )- The effective 
2D interaction strength g characterizes the two-body in- 
teraction, but also depends proportionally on the total 
particle number N and inversely on the oscillator length. 



II. INSIGHTS FROM VARIATIONAL 
CALCULATION 

In this Section we will draw physical insights from 
a time-dependent variational calculation. This is not 
expected to give quantitative predictions, but will pro- 
vide valuable intuition about the stationary solutions and 
characteristic trajectories of the vortex dipole. We use 
the following class of simple variational wavefunctions: 

V = A(t) [z- Zl (t)]{z* -z* 2 (t)] M\z\ 2 ) , (2) 



FIG. 2: (Color online.) Half-distance between pair for sta- 
tionary vortex dipole configuration, function of g. Solid 
line: from Thomas-Fermi ansatz; dashed line from Gaussian 
ansatz; filled circles from time-independent GPE. The dot- 
ted line shows the Thomas-Fermi boundary of condensate. 
Lengths are plotted in trap units. 

where z = x + iy is the complex 2D position coordinate, 
z\ and Z2 are vortex and antivortex positions, and A(t) is 
a normalization constant. Both Gaussian and Thomas- 
Fermi forms are used for the condensate shape function 
/ c . Details of the formulation are in the Appendix. 

Our form ([2]) is amenable to exact treatment. Its dis- 
advantage is that the vortices are 'too big', being deter- 
mined by the trap size rather than the healing length 
£, which is reasonable only for small interactions. We 
note however that in two dimensions £ ~ g~ x ^ decreases 
rather slowly with the interaction parameter, so that vor- 
tices remain big even for relatively large values of g. For 
the physics of stationary states, our wavefunctions there- 
fore give good answers up to large g (up to g ~ 90). The 
situation with dynamics is more complicated. 



A. Stationary solution 

We restrict ourselves to initial positions of the vortex 
and antivortex which are symmetric around the trap cen- 
ter; without loss of generality we place both on the 2-axis, 
yi (0) = y 2 (0) = 0. If we start with x 2 (0) = -xi(0) < 0, 
mutually-driven (inhomogeneity-driven) motion would 
tend to drive both defects in the negative (positive) y 
direction. Mutually-driven (inhomogencity-driven) mo- 
tion wins when the vortex and antivortex are close to 
(far from) each other. Fig. [T] displays this competition 
by plotting the initial velocity, yi(0) = y^iO), against the 
initial dipole separation, x^ — 2xi(0). The initial veloc- 
ity is negative (positive) for smaller (larger) x^. 

We also note that the vortex dipole initial velocity ap- 
proaches that of a single vortex, placed at xi(0) — Xd/2, 
when the distance is large enough, demonstrating that 
the motion of the two defects are indeed determined pri- 
marily by the trap shape when they are far enough apart. 

At some intermediate distance, the two effects balance 
each other, y' x 2 (0) = 0. This nontrivial stationary config- 
uration is a central theme of this work. We denote as x s 
the stationary value of x%(0) = x&/2. We explain later 
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FIG. 3: (Color online.) Vortex and antivortex trajectories, 
calculated using the variational formalism with wavefunction 
(|A.3|) . The four pictures show trajectories resulting from ini- 
tial positions xi(0) = 0.3, 0.6, 0.9, 1.2. Full and dashed lines 
show g = 6 and g — 40 respectively. Lengths are plotted in 
trap units. 



that above a critical interaction g — g c ~ 18.0, a vor- 
tex dipole configuration is a true stationary solution of 
the full Gross-Pitaevskii solution such that not only the 
singularity positions but also the density is stationary. 

Fig- 121 shows the stationary point positions func- 
tions of the interaction strength <?, calculated using the 
trial wavefunctions (|A.2[) and (|A.3|) . Also shown are nu- 
merically exact values of x s , calculated directly from the 
GPE fSec. IIV[) for g > g c as 18.0. The comparison shows 
that the gaussian ansatz is not useful in the relevant re- 
gion g > g c . The Thomas-Fermi ansatz (|A.3|) gives good 
results for the static solutions until g ~ 90. We will there- 
fore restrict ourselves to this wavefunction in Sec. Ill Bl 

For very large g, the ansatz (|A.3|) yields decreasing 
behavior for the x s (g) function, presumably due to the 
'large- vortex' nature of our wavefunction. 



B. Trajectories 

The full curves of Fig. [3] are vortex pair trajectories, 
calculated using the Thomas-Fermi wavefunction (j A.3|) 
for g = 6. When placed close together (cci(0) < x s ), the 
pair starts by moving in the negative y direction, due to 
mutually-driven dipole motion. As they move away from 
the x-axis, inhomogeneity-driven motion becomes more 
and more important, and eventually they turn away from 
each other and move in the positive y direction, crossing 
the x-axis at some x\ > x s . As they move in the +y 
direction, the distance between them decrease, and this 
causes the direction to reverse again. The reversals of 
motion might be regarded as 'reflections' of the propa- 
gating vortex dipole, at each 'end' of the trap. 

The variational calculation gives periodic orbits for 
vortex and antivortex (first three panels of Fig. [31 left 
panel of Fig. 0|. When the starting positions x%(0) are 
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FIG. 4: (Color online.) Vortex position (xi,yi) as function of 
time. Left panel: variational calculation with Thomas-Fermi 
ansatz (g — 6). Right panel: non-interacting bosons solved 
in terms of harmonic oscillator wavefunctions. All quantities 
are plotted in trap units. 



moderately larger than x s , a similar closed orbit pair re- 
sults, starting from the outer edges of the closed trajec- 
tories. The closed orbits are smaller if the initial distance 
from the stationary point, |xi(0) — x s |, is smaller. The 
dashed lines and the last panel in Fig. [3] show unphysical 
features in trajectories with large g or large £i(0), which 
are artifacts of our simple variational setup. 

We will show in Sec. [V] that the full Gross-Pitaevskii 
solution supports very similar trajectories for large g. 
Our simple variational approach has thus already uncov- 
ered a very characteristic motion of vortex dipoles in 2D 
trapped condensates, namely, trajectories in which the 
two defects each revolve around a stationary point. 



III. NON-INTERACTING BOSONS 

The case g — 0, i.e., non-interacting bosons in a har- 
monic trap, is exactly solvable. One can build any desired 
initial state as a linear combination of appropriate har- 
monic oscillator eigenstates. The evolution of each eigen- 
state is known — an eigenstate with energy E evolves as 
4>(t) = <fi(0)e~ lEt . Therefore the evolution of any state 
can be found once it is expressed in this basis. 

A state will be stationary only if it is built out of eigen- 
states with the same energy E. Considering the wave- 
functions of a 2D harmonic oscillator (e.g., Ref. [2fj|). 
it is easy to convince oneself that an off-center vortex- 
antivortex pair cannot be formed out of the eigenstates 
corresponding to a single energy. Thus a stationary vor- 
tex dipole configuration does not exist for g = 0. 

However, one can form a quasi-stationary state in 
which the defect positions are stationary, although the 
condensate density undergoes complicated periodic mo- 
tions. This is formed by combining eigenstates from the 
three lowest eigenvalues, so that there are off-center vor- 
tex and antivortex at (±xi(0),0). It is straightforward 
to show yi(0) = xi(0) - l/xi(0), also plotted in Fig. Q] 
The quasi-stationary configuration is therefore x s = 1 . 

The defect trajectories for g = 0, found from the 
above linear combination, involve sharp direction rever- 
sals along the same path. This suggests that the vortex 
picture is not very suitable here. In fact, the density pro- 
files do not show well-defined vortex and antivortex, but 
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FIG. 5: (Color online.) The stationary soliton solution bifur- 
cates at g«18, producing a stationary vortex dipole solution. 
Dashed (full) line shows energy of a stationary soliton (sta- 
tionary vortex dipole). 



a 'soliton'-like band of low density, even when the system 
contains a dipole-like pair of singularities at (±xi(0),0). 

The right panel of Fig. [4] shows the defect trajectories 
for the 2D non-interacting trapped condensate, via the 
functions x\(t) and yi(t). The x\{t) curve has half the 
period of y\ (t) , indicating that the defects move back and 
forth along the same line instead of around a closed orbit. 



IV. STATIC SOLUTIONS: SOLITONS & 
VORTEX DIPOLES 

The absence of a vortex-antivortex stationary solution 
for g = 0, and some of the dynamics to be presented in 
Sec. |V1 can be put into context by considering a bifurca- 
tion phenomenon (l3l | . 

At small <?, we do have a stationary solution: a dark 
soliton-like object instead of a vortex dipole. This is easy 
to see for g = 0: the first-excited (n x ,n y ) =(0,1) eigen- 
state, Voi oc ye~^ x2+y2 ">/ 2 , can be regarded as a dark 
soliton analog. The term 'soliton' comes from the anal- 
ogy to the case where the confinement in the y direction 
is much weaker, i.e, a 'quasi-lD' condensate. Although 
the word evokes a picture of a propagating object, in a 
trapped condensate a soliton can be stationary. 

As g is increased, the soliton-like configuration remains 
a stable stationary object until g — g c , at which point it 
bifurcates into two solutions, of which the unstable one 
is still soliton-like and the stable one is the vortex dipole 
configuration. Using a many-variable Newton-Ralphson 
method to find solutions of the time-independent GPE, 
we have found that the point above which the stationary 
vortex dipole exists is g c w 18. This is consistent with 
Ref. [13J, but provides better precision. The energies of 
the soliton solution and the vortex dipole solution are 
shown in Fig. [5] as a function of g, and density profiles 
for stable stationary solutions are in Fig. [6l 

The same bifurcation happens in elongated traps, 
where the word 'soliton' is more appropriate. In a 3D 



FIG. 6: Density profiles of stationary states: (1) soliton at 
g = 11; (2) vortex dipole at g — 25; (3) vortex dipole at 
5 = 60. 

trap one gets via bifurcation a stationary vortex ring in- 
stead of a vortex dipole Q . 

The filled circles in Fig. [5] indicate stationary vortex 
dipole positions for g > g c . Note that the value of x s 
is near the trap oscillator length (x a ~ 1) over a wide 
range of g. This is well below the Thomas-Fermi size 
of the system (dotted line in Fig. [2J, which increases as 
i?TF ~ ff 1 ' 4 with the interaction. 



V. NUMERICAL EVOLUTION OF 
GROSS-PITAEVSKII EQUATION 

We now present results from explicit numerical evo- 
lutions of the 2D time-dependent Gross-Pitaevskii equa- 
tion. In these simulations, a vortex dipole was first cre- 
ated at x = ±Xd/2 by propagating in imaginary time, 
at each step normalizing and also imposing the phase 
structure of (z — z\){z* — z%)- After reaching a stable 
dipole configuration, propagation in real time was fol- 
lowed. Time propagation was performed using the 4th- 
order Runge-Kutta algorithm, with the kinetic energy 
term calculated in momentum space. 

Since we have no dissipation, the energy should remain 
constant. Breakdown of the numerical scheme, which 
always happens after some amount of propagation time 
due to finite spatial and temporal grids, is signaled by 
a rapid increase of the energy. The simulation breaks 
down earlier for larger values of interaction g, because 
the length scale £ ~ <? -1 / 4 is smaller. With spatial grids 
of up to ~ 100 per trap oscillator length, we obtained 
useful results up to g < 200. 

Locating the vortex and antivortex is tricky, especially 
for smaller g. It is then often not sufficient to analyze 
the densities. One possibility is to observe both real and 
imaginary parts of ip and find points where they both 
vanish; another tactic is to track phases and identify pla- 
quettes around which there is a 2ir phase winding. 

Our main observation is that the motion of the vortex 
dipole is never periodic, as opposed to the variational 
description. For small <?, the motion is far more compli- 
cated. For larger g, the vortex pair starts to make almost 
closed orbits. Fig. [6] shows one of these relatively easier 
cases, at g — 150, where we have followed the dynam- 
ics for slightly longer than one 'period'. As in the vari- 
ational calculation, the competition between mutually 
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FIG. 7: (Color online.) Vortex pair trajectories from nu- 
merical solution of time-dependent GPE, g = 150. Starting 
positions are a;i(0) = 0.8, 1.1, 1.3, and 1.5. Successive vor- 
tex/antivortex positions have been joined by straight lines. 
The condensate radius (Thomas-Fermi) is _Rtf — 3.72, not 
visible in these plots. Lengths are plotted in trap units. 
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FIG. 8: (Color online.) Vortex pair trajectories from numer- 
ical evolution of time-dependent GPE, g = 10 and g — 50. 
Starting positions are £1,2(0) = ±0.8. Successive positions 
have been joined by straight lines. There are some gaps in 
time (widely separated points joined by a straight line). These 
correspond to intervals where the presence of several vortex- 
antivortex pairs made it difficult to meaningfully track a single 
pair. The Thomas-Fermi value for the condensate radius are 
Rtf — 1-89 and Rtf =— 2.825 in the two cases, shown with 
dotted lines. Lengths are plotted in trap units. 



driven and inhomogeneity-driven motion is clear, as is 
the tendency to reverse direction when the vortex dipole 
reaches the 'edge' of the condensate. There are several 
additional unexplained features. One is the structure at 
the outer edges of the orbits. Also, for large initial pair 
distances iEi(0) > x s , we note that the defects first move 
in the — y direction for a very short time, before starting 
in the +y direction. These two features might be due to 
effects of the relatively nearby boundary [2^, [28| . 

The situation is significantly more complicated for 
smaller g, as shown in Fig. [8] The vortex positions 
are joined by straight lines. The big jumps in these 
'trajectories' represent intervals where the singularities 
were difficult to track or distinguish, or there were ad- 
ditional singularity pairs. The vortex and antivortex in 
the g = 50 case each make almost one complete revolu- 
tion around the stationary point, before coming too close 
to each other, and after this the situation for some time 
becomes too complicated to track. This is when the pair 
undergoes some type of 'annihilation', leading to some 
complicated excitations including additional singularity 
pairs. When we next manage to track a pair, they move 
again roughly along trajectories around the stationary 
point, until they come near each other again. 

The g — 10 case is even more complicated. For g < 18, 
there is no true stationary vortex dipole solution; there- 
fore the paradigm of each defect rotating around a sta- 
tionary position is not expected to provide a useful de- 
scription. Nevertheless, some remnant of the pair 're- 
volving around stationary positions' seems to be visible 
in the defect trajectories. 



VI. DENSITIES 



It is only at rather large <?, when the vortex size £ is 
small, that density plots are useful for identifying vor- 
tices. The middle figure in Fig. [5] shows that for g mod- 
erately above g c « 18, the density profiles for stationary 
vortex dipole configurations look substantially like the 
stationary soliton case (Fig. [6] left). Obviously in a dy- 
namic situation at such interactions, the vortices retain 
little of their individual identity as they move around, 
which helps explain why the singularity trajectories in 
Fig. [8] are not easy to explain in the vortex-antivortex 
language. 

For g~50-100, the vortices are more well-defined in 
stationary density profiles (Fig. [6] right), but still have 
significant shape dynamics when they are moving. It is 
only at very large g that the vortex and antivortex can 
be regarded as point objects. This is the case where the 
Lagrangian analysis, based on wavefunctions with well- 
defined vortices, gives good insight. 

Density plots of the non-interacting Bose condensate 
with a 'vortex dipole' (Sec. typically show a 'soliton'- 
like profile rather than distinct vortex and antivortex, 
indicating that the vortex picture is of limited relevance. 

Our variational wavefunctions (|A.2|) and (|A.3|) also 
have large vortex sizes. As a result, when the pair dis- 
tance is not very large, the density profile tends to show a 
large stripe of low-density region (i.e., soliton-like profile) 
instead of individual depressions for vortex and antivor- 
tex. However, the distinct identities of the vortex and 
antivortex (as well as lack of additional dynamics possi- 
bilities) is built into these wavefunctions. 
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VII. DISCUSSION 

To summarize, we have analyzed, using several 
methods and viewpoints, the trajectories of a vortex- 
antivortex pair in a circularly trapped Bose-Einstein 
condensate, described by the Gross-Pitaevskii equa- 
tion. The principal intuition we have provided is in 
terms of the competition between mutually-driven and 
inhomogeneity-driven motions of the pair. This leads to 
a simple understanding of the stationary vortex dipolc 
solutions [HI, and also leads to the expectation of 
a type of motion where the vortex and antivortex each 
revolve around the stationary positions. Although our 
time-evolution simulations do not extend to interactions 
beyond g = 200, our available results strongly indicate 
that this type of pair motion becomes more and more 
relevant for larger interactions. 

The fact that low-g condensates show more compli- 
cated vortex dynamics has been given proper context 
through our considerations of the non-interacting case, 
the stationary solutions below and above g = g c ~ 18, 
and density profiles, in Sees. IIII1 IIV1 IVII At large g, the 
vortex shape distortion dynamics is less important (be- 
cause they are smaller objects), there is lesser tendency of 
pairs to morph into or mimic soliton-like objects, extra 
pairs are not created as easily, and annihilation is less 
likely in the absence of dissipation. All this favors the 
type of motion described above, obtained from our varia- 
tional wavefunction which only allows position dynamics 
for a single vortex dipolc. At smaller g, the description in 
terms of vortex positions alone is less adequate, because 
of additional dynamics possibilities. 

Vortex dipoles can be regarded as the two-dimensional 
analogs of dark solitons in one dimension, and of vor- 
tex rings in three dimensions. Each of these are soli- 
tary waves in uniform (untrapped) Bose condensates 
of the appropriate dimension, and can become non- 
propagating stationary objects in trapped condensates. 
Vortex rings in trapped 3D BECs display self-propelled 
and inhomogeneity-driven motions in opposite directions, 
just as we have found for vortex dipoles in 2D. 

Our study opens up a number of additional questions. 

1. As described in Sec. El there are aspects of the trajec- 
tories that we do not understand even for relatively large 
<?, such as the initial negative-?/ motion for x±(0) > x s 
(Fig. [7]) . Presumably, this requires understanding of the 
coupling of the vortex position dynamics to the vortex 
shape dynamics. It is intriguing to ask whether the 
trajectory becomes smoother and more periodic at even 
larger g. 

2. For small g, an additional open question is how 
the vortex dipole motion is affected by the presence of 
'nearby' solitonic solutions, e.g., this might explain the 
appearance of extra pairs. (A soliton may be regarded 
as a very large number of vortex dipoles.) 

3. The 'reflection' of the dipole when it reaches the 'edge' 
of the condensate could be better studied in the setting 
of an elongated 2D condensate, where edges are better 



defined than in our circular case. 

4. There remains the open question of vortex dipoles 
which arc initially not placed symmetrically around the 
trap center. 

5. We have not studied the effects of dissipation on vor- 
tex dipole dynamics. In particular, dynamics of vortex 
pair creation and annihilation is important for Kosterlitz- 
Thouless physics [2j. At present we do not have a satis- 
factory theoretical framework for including such effects. 

We hope to address some of these issues in future work. 



APPENDIX: VARIATIONAL FORMALISM 

We describe here the variational Lagrangian formula- 
tion used in Sec. [TT] to follow the motion of vortex and 
antivortex positions. Our results are based on variational 
wavefunctions of the form where vortex positions ap- 
pear as time-dependent variational parameters Z\ (t) and 
Z2(t). Using such a wavefunction in the Lagrangian 



L= dr 



dip dip* 



+ ^*V 2 ^ -V tI (r) H 2 -\g |V| 4 



(A.I) 



one can derive Euler-Lagrange equations of motion for 
the complex parameters z±(t) and z-z(t). Solving the 
equations of motion then provides the trajectories of the 
topological defects. This method for vortex dynamics in 
BECs was pioneered in Ref. m a nd has since been used 
in several vortex applications [2ll. [22l . [23l . [2~H . [25l [2(| . 

For the condensate shape function f c of Eq. @, we 
used both Gaussian and Thomas-Fermi forms: 



^ G = A G (t) [z-zi{t)][z*-zZ(t)] exp[-H 2 /2l 



and 



V'TF = A TF (t) [z-Zi(t)][i 



TF 



(A.2) 



(A.3) 



with /tf = 



a- 5 2 



j g and /x 



Here A(t) is a time-dependent normalization constant. 
Retaining this factor A(t) is essential for getting the dy- 
namics even qualitatively correct. For example, for the 
single- vortex case, omitting this factor gives vortex pre- 
cession in the reverse direction! Since A(t) actually con- 
tains the vortex position parameters zi^it), terms in the 
equations of motion are missed when A(t) is omitted. 

Our form ([2]) is generalized from the lowest-Landau- 
level form for vortices in a rapidly rotating condensate 
f27l ]. It has the advantage of being amenable to exact 
treatment, i.e., the integrations in (|A.1[) can be per- 
formed analytically. The disadvantage is that the vor- 
tices are 'too big', as discussed in Sec. |TTJ We could in 
principle use any wavefunction of the form 



1> = A{t) g v (ui)e 1 ^ SvMe-* f c (\z\ 



(A.4) 



7 



where Ui = \z — and (f>i = tan 1 { v x _ v x i . )■ The vortex 
shape function g v (u) ideally would be linear in u for u — » 
and constant for large it Q; for example 

V" + ? 

Eq. ([2]) corresponds to g v (u) — u, having the wrong 
behavior at large u. This is why our variational wave- 
functions have unnaturally large vortices. Unfortunately, 
we were unable to identify a function g y (u), having the 
correct limiting behaviors, for which the integrals in 
Eq. (jA.ip can be done analytically. We also did not 
attempt to incorporate velocity field effects due to the 
boundary [2^, [28|. Our presented results are therefore 
limited to wavefunctions of the form of Eq. 

Using one of our variational wavefunctions, (|A.2|) 
or (|A.3|) . the Euler-Lagrange equations of motion, 
D t (dL/du) = dL/du, become equations of motion for 



the defect positions Xi^t) and yi^(t)- We can obtain 
these equations of motion analytically, but they are too 
cumbersome to reproduce here. 

For a single vortex-antivortex pair with unrestricted 
positions, there are four variables, leading to four cou- 
pled first-order nonlinear differential equations which can 
be solved numerically for the vortex trajectory. For the 
initial conditions we have used, we expect the motion 
to be symmetric around the y-axis, so that we can use 
x\(t) = —x^it) and yi(t) — J/2 (*) • This significantly re- 
duces the numerical cost of solving the equations of mo- 
tion, of which there are now only two. 
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